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A TECHNIQUE FOR IDENTIFYING PILOT DESCRIBING 
FUNCTIONS FROM ROUTINE FLIGHT-TEST RECORDS 
By Rodney C. Wingrove and Frederick G. Edwards 
Ames Research Center 


SUMMARY 


Previous studies have shown that the dynamic response of the pilot can be 
represented by a linear element (describing function) and a remnant term 
(output noise) . The previous work also has indicated that there is an error 
in identifying the pilot describing function from routine tracking task 
records because the output noise of the pilot transfers through the control 
loop, producing an undesired correlation with his input signal. This report 
shows that this correlation, and thus the identification error, can be 
reduced by shifting the input signal during the computer processing an amount 
equivalent to the effective time delay of the pilot. This report includes a 
theoretical analysis of this technique and examples to illustrate its 
application. 

The theoretical analysis considers the fact that the computer processing 
is constrained to identify only physically realizable systems. With this con- 
straint, it is shown that the error in identifying the pilot describing func- 
tion depends on the spectrum of the pilot's output noise; the identification 
error can be made small if the noise is near "white" in relation to the sum 
of all effective time delays through the control loop (pilot plus controlled 
element). This result is significant because, if these conditions are met, 
it is possible to identify the describing function of the pilot in a feedback 
system that is excited only by his output noise. 

The identification of several simulated pilot models is included in this 
study to illustrate this technique. Also, representative data from the retro- 
fire phase of the Gemini X flight have been analyzed and are presented to 
demonstrate the successful application of this technique with routine 
spacecraft operating records. 


INTRODUCTION 


This report considers the problem of identifying the input-output rela- 
tionship of the pilot by use of measured data from routine flight operations 
in which the pilot provides feedback control. The problem in using the mea- 
sured input and output data directly is that any extraneous output noise by 
the pilot causes an error in identification. This problem is solved in this 
report by the development of a computer processing technique that, under cer- 
tain conditions, yields an estimate relatively free from identification error. 



Since the identification of feedback control systems is important in many 
fields, the technique has wide significance and applicability. 

The input-output characteristic of a pilot must be regarded as random, 
nonlinear, and dependent on the task he is to perform. Many previous studies 
have shown that this type of response can be represented appropriately with a 
quasi- linear system modeled by a linear element (describing function) and a 
remnant term (output noise). The pilot describing functions usually have been 
identified from records obtained in ground-based simulators (ref. 1) and 
flight tests (ref. 2) wherein carefully controlled external forcing functions 
are used to excite the pilot-vehicle system. The pilot describing functions 
are measured by comparing the input and output signals of the pilot with the 
known forcing function. This method minimizes those errors in identification 
due to any correlation of the input signal with the pilot's output noise. 
Reference 3 contains a good review of this previous work and summarizes the 
measured pilot describing functions. 

Most other methods for measuring pilot describing functions depend on 
random disturbances (e.g., aerodynamic turbulence, propulsive disturbance, 
etc.) to excite the pilot-vehicle system. These methods compute directly the 
describing function of the pilot from his input and output signals. However, 
there is a fundamental difficulty with these methods because the pilot's out- 
put noise transmitted through the control loop produces an undesired correla- 
tion between his input and output signals, thereby causing an error in 
identification. In reference 4, the expected error is analyzed and it was 
shown that if the amplitude of the pilot's noise is large, as compared with 
the external disturbance, then the identification error is unacceptable. 

During routine flight-test operations, there are no carefully controlled 
forcing functions and even the random external disturbance may be quite small 
so that the principal system excitation may come from the pilot's output 
noise. This report shows that in such situations it may still be possible, 
under certain reasonable conditions, to determine the pilot describing func- 
tion without incurring an unacceptable identification error. One required 
condition is that the pilot (or possibly the feedback control loop) have a 
time delay. If this condition is met, it is possible to take advantage of 
this fact in the identification data processing. In effect, the input signal 
is shifted during processing by an amount equal to the time delay of the 
pilot. Although previous studies (refs. 5-7) have considered the use of a 
time shift in the measurement of pilot describing functions, it was apparently 
not observed that this time shift would strongly influence the error in 
identification. 

This report presents a theoretical analysis to show that this technique 
will reduce the identification error. The simulation and identification of 
several known system elements are included to compare with the theory and to 
illustrate the use of this technique. Also, results obtained from the retro- 
fire phase of the Gemini X mission are presented to demonstrate the 
application of this technique to routine flight-test records. 
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controller deflection (output of pilot) 
error signal (input to pilot) 

Fourier transform of [ ] 
external disturbance 
constant gain 

numerator terms in Y c (joj) 
internal noise (pilot remnant) 
cross-correlation function of e(t) and c(t) 
cross-correlation function of e(t) and n(t) 
autocorrelation function of e(t) 
autocorrelation function of n(t) 
time, sec 

controlled element describing function 
measured describing function (ideal) 
measured describing function (actual) 
pilot describing function 
estimated pilot describing function 
exponential decay factor, sec " 1 
residual 

time shift used during analysis, sec 

time delay in Y c (joj), sec 

time delay in Yp(joi), sec 

power spectrum of e(t) 

power spectrum of n(t) 

cross-power spectrum of e(t) and c(t) 



cross-power spectrum of e(t) and n(t) 
co frequency, rad/sec 


BACKGROUND 


This section discusses the piloted control system elements and indicates 
the error in identifying the pilot describing function from routine tracking 
task records. A computing process for reducing this identification error is 
then outlined. This background material precedes a more detailed analysis of 
the identification error presented later in the report. 


General Remarks 

Figure 1 is a block diagram of the pilot in a compensatory tracking task 
trying to control his output c(t) so that the input error signal e(t) 


Pilot Controlled system 



Figure 1.- Identification using standard measurement methods. 


is kept near zero. Generally, the input-output characteristics of the pilot 
must be considered as complex, nonlinear, and time varying. However, for the 
purposes of modeling, it is common practice to assume that his characteristics 
can be represented by a quasi-linear system (ref- 3). This mathematical model 
contains the linear element Yp and the noise source n. The element Yp(jw), 
which is called the pilot describing function, 1 is a linear constant- 
coefficient system with a frequency response dependent on the in put e (t ) . 

1 Technically, Yp(jco) represents a random input describing function 
because random^ rather than sinusoidal, signals are used here (see ref. 3). 
Also, to avoid additional notation, terms such as Y(jm) will be used to 
represent both the transfer functions of linear systems and the describing 
functions of nonlinear systems. 
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The term nft) represents the difference between output of the pilot, eft), 
and output of the describing function Ypfjoo) driven by eft). Thus nft) 
accounts for remnant terms such as nonlinearities, time variations, and 
additive noise in the output of the pilot. 

The controlled system is mathematically characterized by the constant 
linear element Y c and the noise source i. The time history ift) accounts 
for nonlinearities and time variations in the controlled element, time- 
varying commands, and all disturbances from aerodynamics, propulsion, etc., 
external to the pilot. 


Identification Error Using Standard Methods 

Several methods (e.g., refs. 4-10) have been used to compute, from given 
records of eft) and eft), a describing function Y m (jeo) that represents the 
best linear relationship between eft) and eft). Best here means that the 
integral of the squared residual, /e 2 ft)dt, is minimized over a given record 
length, where eft) is the difference between the actual record cift) and the 
output of the system Y m fjco) excited by eft). The measurements Y m (jco) may 
differ somewhat between methods because each method uses slightly different 
approximations and model forms in computer processing. Generally, the measure- 
ments of V fjco) will be near the following ideal describing function Y m (jD 
that represents the best linear relationship between eft) and eft) for random 
stationary signals 2 : 


Y m (jo0 « Y m (jw) = - CD 

$ ee CD 

In this equation, $ ec fjw) is the cross-power spectrum between eft) and 
eft) and $ ee (co) is the power density spectrum of eft). In identifying the 
pilot describing function with these types of methods, previous studies (e.g., 
refs. 4 and 9) have shown that there is a difference between the measured 
describing function Y m (ja)) and the actual describing function Ypfjco). This 
difference, or "identification error," can be shown by delineating the compo- 
nents of the cross-power spectrum: ^ecU^) = Yp (jeo) $ ee (u) + $en(iu). Sub- 

stituting these components into equation (1) yields 

$ en (j D 

v n (j») - YpCM * (2) 

error 

Equation f2) shows that any cross-correlation $en(j w ) will contribute an 
error in identification. Such a correlation does exist during closed- loop 
control because nft) transfers through Y c fjco) and thus appears as a compo- 
nent of eft). If nft) is much smaller than ift), the ratio $en (jw)/$ ee (co) 
will be small and the measured transfer function Y m fjco) will be near the true 


2 If the measurement has the constraint to identify only physically real- 
izable systems, then, as shall be pointed out later, Y m (joo) is written in a 
slightly different form. 
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value Yp(jw). However, if n(t) is much larger than i(t), the ratio 
$en(J“)/ $ eeC u ) will be significant and the measured describing function 
Y mO w ) wil l be very different from Yp(ju)). For routine flight-test condi- 
tions, where n(t) may be much larger than ift), it is necessary to find some 
means of reducing this error. Such a technique will be outlined next. 


Use of a Time Shift in Identification 

Previous studies (e.g., refs. 5-7) have considered the use of a time 
shift during the computer processing to account for the effective time delay 
of the pilot. This time shifting represents only a slight modification to the 
identification methods in figure 1. 

This time-shifting technique illustrated in figure 2 involves the follow- 
ing steps in the computing process. 

1. The input signal e(t) is shifted with respect to c(t) by an amount 

A, where A is equivalent to the time delay of the pilot. 

2. The describing function Y m (ju) is determined using the shifted data 

from step (1). 

3. The estimated transfer function is determined from the measured 
transfer function as 


YpCjoO = e" Xj “y m (joO 


Pilot 


Controlled system 



Figure 2.- The use of a time shift X in identification. 


Although previous studies have considered this time-shifting technique, 
it was apparently not observed that this technique would strongly influence 
the errors in identification. This report shows that when this technique is 
used with a measurement method in which is constrained to be 
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physically realizable, 3 then the identification error due to the correlation 
of e(t) with n(t) can be reduced. This reduction will be shown in the next 
section where the identification error to be expected with this computing 
process will be analyzed. 


ANALYSIS OF IDENTIFICATION ERROR 


The reduction of the identification error by the foregoing computer 
processing will be illustrated from two points of view. First, a general 
analysis will show why the time shift X reduces the identification error. 

The second analysis will develop equations to show, in more detail, the amount 
the error is reduced. The following analysis is presented using the frequency 
domain. A similar analysis is presented in appendix A using the time domain. 


General Analysis 

To illustrate the reduction in identification error, equation (1) is 
rewritten as 


Y m O) « Y m (jaO 


F [R e c( T )] J.^ec ^) 6 ^ dx 

F[Ree(x)] J*“ R ee (T)e- T 3“ dt 


( 3 ) 



Figure 3.- Effect of time shift on correlation 
functions . 


where F[R 6 C (t)] represents the 
Fourier transform of the cross- 
correlation function R ec (x) anc * 
F[R ee C T )l represents the Fourier 
transform of the autocorrelation 
function R ee (T). Representative 
curves 4 of the measured quantities 
R ec (x) and ReeM are sketched in 
figure 3(a). The error contribution 
R en ( T )> contained in R e c( T )> is also 
shown for comparison. 

Now consider those measurement 
methods that have the constraint of 
physical realizability. These 
methods used only data for positive 
values of x, and, accordingly, the 
measured transfer function is 


3 This constraint is inherent in the computer processing for most time- 
domain measurement methods such as cross-correlation (refs. 4, 8, and 9), 
orthogonal filters (refs. 4, 5, and 7), and parameter trackers (refs. 4, 6, 
and 10). Most frequency domain measuring methods using cross-spectral com- 
puting programs (ref. 1) usually do not contain this constraint. However, 
such a constraint could probably be incorporated. 

4 These data are from example 1, which will be discussed later. 
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Y m<» 


Yp(jt*0 



error 


(5) 


With this constraint, only that portion of R en (x) f° r positive values of x 
(shown by the shaded region in fig. 3(a)) contributes an error in 
identification. 


Let us next introduce the time shift X as presented in figure 2. This 
time shift is applied so that the shifted input data are e'(t) = e(t - X). 
The effect of this time shift is illustrated in figure 3(b) where the func- 



(b) With time shift X 
Figure 3.- Concluded. 


tions R e f nO)> Re'cCO* and R e ' e ' CO 
resulting from the shifted input 
data are presented. It is shown 
that the addition of this time shift 
causes the quantities Re ! cCO> and 
R e f nCO to be shifted by the amount 
X with respect to R e » e t(T). Now it 
is apparent that the error contribu- 
tion of R e t n (r), for the positive 
values of t, is reduced and 
includes only the small shaded area 
in figure 3(b). The actual value 
for the error term is 


J Q R e 'e’ COs"^ dr 


J/enCx + X) e - T ^ dx 
J”Ree(x)e“ T;ia) dx 


( 6 ) 


In general, R en (x) w iH decrease for positive values of t. Therefore, note 

that the error contribution f R en (r + X)e dr will be reduced as X is 

0 

increased. Further, the error contribution will be zero if RenCO is zero 
for values of t greater than X. 

This general discussion has attempted to give some physical insight into 
why the time shift X reduces the error in identification. We will now turn 
our attention to a more detailed analysis to determine the amount the error 
can be reduced. 
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Detailed Analysis 


In this section, we will derive formulas that show the reduction in 
identification error as a function of the primary variables within the control 
loop. As noted earlier, we will consider the use of a time shift X in the 
identification and consider measurement methods in which Y m (jw) is 
constrained to be physically realizable. 

To mathematically represent a measured describing function Y m (jaO that 
is constrained to be physically realizable, we can utilize the relationship 
used with the Wiener-Hopf equation (ref. 11). Using this relationship for 
physically realizable systems, equation (1) is written 


Y m(W 


Y m (ja>) = — 

S’eeCjw) 


4>ecCj“) 

SeeCjw) 


(7) 


where 


$ ee o) = 

^eeCjm) 

[ U 


i’eeOwD^eeCjw) 

has poles or zeros only in the left-half plane 
has poles or zeros only in the right-half plane 
has poles only in the left-half plane 


This follows the usual form, which implies that the direct transform of a time 
function that is stable and zero for negative time will have all its poles in 
the left-half plane (LHP) . 


Now we introduce the time shift A as illustrated in figure 2 and define 

the shifted data as e T (t) = e(t - X). Because ^e'cU^) = e ^ W$ ec (j w ) and 
$ e » e t(w) = $ ee (u)), we can write the measured transfer function as 


Y mCjw> = 




W $ec 0“) 

. ®ee Ww) - 


(8) 


As shown in figure 2, we can define the estimated describing function in 
terms of the measured describing function, YpCjw) = w ) . And if we 

assume that there is no modeling error, that is, YmCjw) = Y m (jio), then a 
theoretical expression for the estimated describing function is 


Yp(jco) 


e~ x ^ 

$ eeO) 


"e'^ aj $ ec (ja>) 


®ee 0 


(9) 
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Introducing the individual terms for $ ec Cjw) (see eq. (2)), we have 


Y p CjoO 


$ eeCi w ) 


p (jaO$g e (juO 

L r J + 


- X J CO 


^eeCju) 


g Xj^ e n 

$ ee^ u ) 


( 10 ) 


The impulse response function of Yp is assumed to be zero for time less 
than a value of Xp and, so long as \ is less than or equal to Xp, the 

term e^^YpCjoi) has poles only in the LHP. Simplifying equation (10) with 
this assumption, we obtain 


„ - Xj 0) 

Yp(ju)) = YpCju) + -| 

$eeCju) 


e^jw$en Cjw) 




-> + 


( 11 ) 


The term $ ee (io) consists of contributions from two sources: i(t) and n(t). 

The maximum error can be determined by assuming i (t) = 0 (ref. 4). With this 
assumption and using basic closed-loop relationships (ref. 11), let us define 


$L(j“) = 


-Y c (jw)$ nn (jw) 


ee 


1 + Y p Y c (ja>) 




-Y c (- Jtu) <D nn (jw) 

1 + Y Y (-j w ) 
p c J 


®enO“) = 


Yq ( j w) ^nn (i ^nn (l^^ 
1 + Y p Y c (-j(o) 


( 12 ) 


(13) 


(14) 


These definitions assume that Y c (jto) is minimum phase (i.e., contains no time 
delay or zeros in the RHP). The case in which Y c (jco) is a nonminimum phase 
will be illustrated at the end of this section. 


Y c Minimum phase.- From the foregoing assumptions, which cover a variety 
of piloted control situations, we find that 

[e^’^nnCj^U 


Yp(jw) = Y p (j to) - 


— — + Y ( j to) 

e ^ u# nn(» L Y c (jw) P 


(15) 
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where now the error is conveniently expressed as a function of S’nnO 61 )’ the 
excitation noise source. In equation (15), the terms [e^ u $nn(j<<0] + and 
$nn(jio) can be evaluated as shown in the following equation: 


YfjuO = Y Hco) - 


uu 

J> 


0 nn 


(t + A)e dx 


f 


WOO® t;,oj d- r 


Y c (jw) 


+ Yp(jw) 


~ Aj Co 


(15a) 


The contribution to the error term includes that portion of R n n( T ) f° r va l ues 
greater than X (shaded area in fig. 4). It is seen that this contribution to 

the error term will be reduced as X 
is increased. (However, this theo- 
retical derivation holds only for 
those values of X less than or 
equal to Tp.) 

Equation (15a) indicates those 
conditions under which the identifi- 
cation error will be small. For 
instance, note that if X is posi- 
tive and if n(t) is white noise 

_ ( R nn( T ) is an impulse at x = 0) , then 

I«nn(T + X)e dx is zero and there will be no error in identification. 

More generally, the identification error will be zero if 



Figure 4.- 


Reduction of error contribution with 
A . 


^nn 


where 


for x > A 


(16) 


- T P 


This result appears to have significance for many applications. The most 
important point is that when these conditions are met, a describing function 
Yp(ju)) within a feedback system can theoretically be measured with the system 
excited only by the internal noise n(t). 


In most realistic situations, R nn (x) will not be identically zero for 
values of x > A. We will next show, however, that the identification error 
can be reduced, and in some cases be made quite small, with more realistic 
forms of RnnW- F° r example, assume that the noise n(t) takes the form 

R nn ( T ) = Ke a i T ! which, for small a, would be narrow-band (nonwhite) noise. 
This form agrees quite well with some experimental measurements of the pilot 
remnant. (For instance, this exponential form with a = 5 sec agrees with 
the measured n(t) in references 3 and 12.) With this form, we can evaluate 
the constant factor of equation (15) as 


fV Cw) V j » T dx 






-ax -jeon: 


= e 


-aX 


dx 


and arrive at 
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Y p O) = Y p (j<u) - e 


(17) 


-aX 


_Y c (ju) 


+ YpOco) 


-Ajio 


The error term on the right side of the equation is a function of the magni- 

ct A 

tude of the constant factor e . As \ increases and if a is large (near 
white noise), then Yp(jw) =“Yp(juj). Conversely, if X = 0, then the result 
is identical to that shown in reference 4: Yp(jcu) = -1/Y c (jw). 

Y c Nonminimum phase.- Let us define the nonminimum phase terms as 

e to represent any pure time delay in Y c (jaj) and N~(j w) to represent 

any RHP zeros in Y c (jw). Then, by including these terms, in equations (12) to 
(14), the estimated describing function becomes 


YpO) = Yp(ju) 


le (A+Tc) j ( j <*») N~ (- j co) /N“ ( j co) ] + 

e CX+Tc) j V_ ( j «o)N" (- j co) /N ' (j <o) 


1 

Y c (jco) + 


Y p(i“) 


(18) 


Note that, in this case, if R nn (x) = 0 for x > x C5 then Yp(jw) need 
not have a time delay (and no time delay, x, is required in the analysis) in 
order that the identification error be zero. 

In this more general case, the identification error in equation (18) will 
be zero if 

R nn^ T ^ = 0 for T > x + T c ( 19 ) 


where 


A i t p 


The identification error can be made small if the autocorrelation function of 
the internal noise R nn C^) negligible for values of x greater than the sum 

of the time delays X + x c . We can also note that the term N c (- jaO/N (jw) is 
representative of a Pade" approximation to a time delay. Thus, any RHP zeros 
in Y c (joj) will tend to act as an additional effective time delay and will 
further reduce the identification error. 

This analysis of the identification error indicates that the internal 
noise n(t) need not be a hindrance to identification, but rather it will aid 
in the identification of feedback control systems if the conditions of equa- 
tion (19) are met. This analysis also may have application in many other 
fields such as biology, economics, and chemical processes. Although these 
other applications are not considered in this report, they do contain time 
delays and some of the measurements can be made only with the noise introduced 
within these systems to be identified. 
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APPLICATIONS AND DISCUSSION 


The use of the computing technique outlined in this report will be 
illustrated through the identification of two examples using simulation data 
and one example using actual flight data. Each example will illustrate a dif- 
ferent point. With example 1, the foregoing theoretical results will be com- 
pared with experimental results. With example 2, a method for selecting the 
time shift X will be illustrated. With example 3, an application using 
actual flight records from Gemini X will be illustrated. 

The simulated systems for the first two-examples are shown in figure 5. 
The dynamics for these examples were simulated on a digital computer. The 
output of a random noise program was appropriately filtered to obtain the 
desired spectrum of n(t). The resulting dynamic records of e(t) and c(t) 
were processed using the method described in appendix A. The experimentally 
determined Yp(jto) to be presented for these simulated examples represents the 
average values obtained from 12 separate 20-second runs. 



(a) Example I 



(b) Example 2 


Figure 5.- System examples used to illustrate identification technique. 
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Example 1: Comparison of Theory With Experiment 


To illustrate the theory, the system in figure 5(a) was simulated and an 
identification was made on the known model. The pilot model and controlled 

element were Y p (jw) = 4^~°* 3:]a) and Y c (joo) = 1/jio. The measurements were made 

with no external disturbance, i(t) = 0, and with the only excitation being the 
internal noise source, n(t). Two forms of the noise spectrum were considered: 
an n(t) with a spectrum that approximates white noise to illustrate the con- 
dition from equation (16) for no identification error, and with a spectrum 

whose autocorrelation function is Rnn^ T ^ = e to illustrate the theoret- 

ical results from equation (17) for an expected identification error. 


Internal white noise .- For this case, the excitation source n(t) had a 
spectrum near white noise. The time shift used in the computer processing was 
taken at X = 0.2 sec. These conditions meet those specified for equation (16) 
According to equation (16), the estimation technique will identify the actual 

system, Y p (jio) = 4e"°' 3; j“. 

Figure 6(a) presents the experimentally determined magnitude 
and phase angle i Y (joj) as functions of frequency. Also shown for 


20 - 


5 15 - 


Y (j«*d 
P 

Y p (i“) , 


Actual 


Experiment 


(0 - 


£ 

<r 



{ a ) n (t ) ~ white noise 


Figure 6.- Identification of example 1; 
X = 0.2 sec. 
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Y p (jm)| 
comparison 
are the magnitude |Yp(jaj)| and phase 
angle {Yp(jaj) of the actual system. 
The estimated amplitude |Yp(jaj)| 
varies ±0.5 dB about the actual 
value for frequencies to about 
9 rad/sec and the phase angle <fYp(joO 
is within ±0.5° of the actual value. 
These differences appear to be within 
the experimental accuracies of the 
simulation. These results substan- 
tiate the theoretical conclusion that 
it is possible to identify the 
describing function of a system that 
is excited by noise n(t) introduced 
within the system. 

Internal nonwhite noise .- This 
case uses the same control elements 
and the same value \ = 0.2 sec as in 
the previous case. However, the 
assumed noise spectrum has a more 

'5 I T I 


realistic form R nn (r) = 

this case, the theory (eq. (17)) 
predicts the following estimated 
describing function : 


For 



YpO) 


-0 . 3ju)-^-0 . 2Jaj 
y 

error 


( 20 ) 
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the 


This 


Y„(j 


theoretical value of Yp(jw) is presented in figure 6(b) along with 


w) 


y tjw) ^ 
Y p lj“> ( 


Theory 


Experiment 


obtained from the experimental data. Also shown for comparison 

are the describing functions of the 
actual system Yp(j to) and the negative 
inverse of the controlled element 
-l/Y c (jto). The experimentally derived 
Yp(jm) in this figure is close to that 
predicted by the theory. We can see 
from this figure that the estimated 
magnitude |Yp(jto)| is about 4 dB below 
the actual value, |Yp(jio)|, at the 

lower frequencies and tends to give 
the appearance of lead 
(slope = 20 dB/decade) at the higher 
frequencies. Overall, the estimated 
magnitude tends toward |l/Y c Cj<o)| as 
predicted by equation (17). The esti- 
mated phase angle, however, agrees 
quite well with the actual value. 




Frequency , uj , rad / sec 

(b) R nn (r) -- <? " 5 l r| 
Figure 6.- Concluded. 


Y li ta d. Actual 
Yp^, X = .3 sec 



(a) Describing functions for X=.05sec and X=.3sec 


Figure 7.- Effect of X on identification; 
example 1; R nn (T) = 6T 5 I t '. 


If a time shift were not used in 
this example, that is, if X = 0, then 
the estimated describing function 
would be Yp(jcu) = -1/Y c (jm), shown by 
the line in figure 6(b). It is inter- 
esting to note that with X = 0, the 
value of the constant factor in the 
error term of equation (17) is 

QJ Tl 

e =1. The value of this constant 
factor with X = 0.2 sec, as shown in 

equation (20), is e ~ 0.37. There- 
fore, for X = 0.2 sec the magnitude 
of the error term was reduced 
approximately 63 percent. 

Figure 7 illustrates the effect 
of X on the identification error. 

The control system and spectrum of the 

excitation noise source, Rn n (T)=e , 

are the same as in the previous case. 
Experimental data, Yp(joo), are shown 

in figure 7(a) for X = 0.05 and 0.3 
sec. Also shown for comparison are 
the actual system Yp(joo) and the 
negative inverse of the controlled 
element -l/Y c (joj). These experimen- 
tal data illustrate the effect of X 
predicted by the theory of equa- 
tion (17). For small values of X, in 
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(b) Identif icotion error for w = I rad/sec 

Figure 7.- Concluded. 


very small at A tsa t . As shown by 
near the time delay F Tp, to minimize 


this case X = 0.05 sec, the estimated 
describing function Yp(jm) tends 
toward -1/Y c (jw). For larger values 
of X, in this case X = 0.3 sec, the 
estimated describing function is 
nearer the actual system describing 
function. 

Experimental data are presented 
in figure 7(b) for several values of 
X. These results are for one value 
of frequency, oj = 1 rad/sec. In this 
figure, the experimental data are com- 
pared with the value for the actual 
system ( | Yp | « 12 dB and i Yp ^-18°). 
This comparison illustrates that the 
identification error decreases, as pre- 
dicted by theory, to a value of 
X = 0.3 sec, which is equal to t . 

(The theory, as noted previously , is 
valid only for X <_ Tp.) The error 

in | Yp | is seen to increase for 
values of X much larger than t . 

This is to be expected because Yp(jco) 

cannot properly model Yp(jco) with 

X > tp . These data^ indicate that the 
minimum error in j Yp | occurs near a 

value of X « Tp . The error in the 

phase angle •( Yp is also seen to be 
is example, the value of X should be 
he error identification. 


We shall now recall an interesting point previously discussed (with 
eq. (2)) that can now be reinforced with experimental data. That is, the mini- 
mization of the identification error, as we have shown above, is not equivalent 
to finding the minimum value for the squared residual e 2 (t) integrated over a 
a given run length. The following table presents the experimentally derived 


Time shift X, 
sec 

Normalized fit term 
/e 2 (t)dt//c 2 (t)dt 

0 

6.02 

.1 

. 19 

.2 

.25 

.3 

.27 

.4 

.28 

.5 

.31 


values of the fit term /e z (t)dt, nor- 
malized with respect to /c 2 (t)dt, as 
a function X. This table shows that 
the minimum value for the fit term is 
not at A = t = 0.3 sec, but rather 
is at X = 0. This is to be expected 
because with X = 0 then 
Yp(jw) = -1/Y c (ja)) and there is essen- 
tially perfect correlation between 
e(t) and c(t). 
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The fact that the minimum value for fit term /e 2 (t)dt appears at X = 0 
can lead to an erroneous interpretation in selecting the best value for X. 

For instance, in previous studies (refs. 5-7), X was selected by using that 
value which gave the best correlation between eft) and c(t) (i.e., the mini- 
mum value for /e 2 (t)dt). This previous method of selecting X is unsatis- 
factory, however, because as we have just noted, if^ n(t) >> i (t) , then the 
best correlation is with X = 0 and, in this case, Yp(joo) = -1/Y c (ju)). An 

alternate method of selecting X will be illustrated by the following example 


Example 2: A Method of Selecting the Time Shift 

The previous example pointed out that the time shift X should be near 
the time delay Tp to minimize the identification error. The time delay Tp 
may be approximately known in some situations (e.g., ref. 3) but, in general, 
its value will be unknown and will depend on the particular piloting task. 

This example illustrates one method of selecting X and will consider 
identification of the system shown in figure 5(b). The excitation noise 
source is the same as used in the previous case, R nn (x) = 0“ 5lTl and i(t) = 0 

however, different forms for the pilot model Yp(jto) = 2 ( j to + l)^" 0 ' 5 ^ and 

l 

controlled element Y r (jco) - — were chosen to demonstrate the iden- 

ju)(ja) + 1) 

tification of more complex dynamics. We will assume in this illustration that 
Yp(jaj) is unknown. The objective will be to estimate the value of Tp and 

then use this value for X to obtain a best estimate for Yp(jai). 

For this example, the following procedure was used to estimate x p and 
thus, select X. 

1. Plot the estimated describing function for a selected value of X. 

2. Determine a transfer function that fits the plot, that is, 

“ j to 

Yp(jco) a* (K x + K 2 jw)e F , etc. 

3. Note the value of estimated Xp from step 2. 

4. Repeat steps 1 through 3 until a value of X approximately equal to 

the estimated Xp is obtained. 

The estimated describing functions for example 2 are presented in 
figure 8 for X = 0.2 and 0.4 sec. Also shown for comparison is the describ- 
ing function for -l/Y c (joj). These experimental data for X = 0.2 and 
0.4 sec follow the trends predicted by theory. The curve for the high value 
of X (as compared to the curve for the lower values of X) tends away from 
“1/Y c (jw). IT was found for these data that any value of X between X ^ 0.3 
sec and 0.8 sec resulted in approximately the same describing function as 

shown for X = 0.4 sec. This estimated describing function can be approxi- 

-T D ja) 

mated by a transfer function of the form Yp(jai) (K x + K 2 ju)e F . The 
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Figure 8.- Identification of example 2. 


estimated Tp values from fitting 
this transfer function to the plots 
(e.g., fig. 8) are presented in fig- 
ure 9 for several values of A. It 
is seen that A is equal to the 
estimated time delay, t , at A « 0.5 

sec. Therefore, X = 0.5 sec should 
be selected for use in this example 
identification analysis. It is seen 
that for this example the method 
works well for estimating the actual 


Figure 10 compares the estimated 
describing function using A = 0.5 
sec with the actual describing func- 
tion for example 2. Both the magni- 
tude and the phase angle of the 
estimated describing function are 
seen to be near these values for the 
actual system. For this case, with 
A = 0.5 sec, the theory of equa- 
tion (17) predicts that the identifi- 
cation error in magnitude will be 

— n \ 

e =0.08. The experimental data 
shown in figure 10 appear to be 
within this 8-percent identification 
error as predicted by the theory. 



Figure 9.- Comparison of estimated x with A; 
example 2. 




Figure 10.- Comparison of Y with Y ; 

example 2. P P 
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Example 3: Flight Test Results From Gemini X 


Flight data from Gemini X were analyzed to illustrate the application of 
this identification technique. In analyzing flight-test data, it is best to 
select a section of the record that contains disturbances external to the 
pilot. As noted earlier in the discussion following equation (2), external 
disturbances will tend to reduce the error in identification. The retrofire 
maneuver is a case in which external disturbances are introduced due to the 
unsymmetric ripple firing of the four retrorockets . The relationship of the 
pilot control task, the jet control system, and the disturbances during retro- 
fire is illustrated schematically in figure 11. 


Pilol 



Figure 11.- Pilot describing function and flight control system; example 3. 

During retrofire, the pilot controls the attitude about each of the three 
axes. There is no control coupling between these axes, and the pilot appears 

Retrofires to treat them as three separate tasks. 



Figure 12.- Time history of yaw control task 
during retrofire. 


Of the three axes, the control about 
the yaw axis contained the best consis- 
tent correlation between attitude 
deviations, e(t), and control stick 
deflections, c(t). A time history of 
the recorded yaw control data is pre- 
sented in figure 12. These control 
data will be used to illustrate the 
measurement of the pilot describing 
function during the retrofire of the 
Gemini spacecraft. (It should be 
emphasized that this was a normal 
retrofire maneuver and the astronaut 
had no prior knowledge of this 
identification attempt.) 

The pilot describing function 
obtained for the data of figure 12 are 
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Frequency, oj , rad /sec 


Figure 13.- Identification of pilot describing 
function; example 3. 


presented in figure 13. Curves of 
magnitude, |Yp(ja))|, and phase angle, 

-(Yp(ja)), are presented as functions 

of frequency for X = 0 and 0.6 sec. 
Also shown for comparison is the 
describing function 5 for -1/Y c (ja)). 
The significance of this line was 
noted previously. The theory predicts 
that for X = 0, the estimated 
describing function Yp(jw) will tend 

toward -l/Y c (jto) as illustrated in 
figure 13. However, for this flight 
situation, Yp(jo)) does not coincide 

exactly with -1 /Y c (jgo) because of 
external disturbances due to the 
firing of the retrorockets and jet 
control system. 

For X = 0.6 sec, the estimated 


describing function Yp(jco) tends 

away from the curve of -1/Y c (jaj). Any value of X from about X = 0.3 to 
0.7 sec resulted in approximately the same describing function as shown for 
X = 0.6 sec. This estimated describing function can be approximated by a 

- T pjw 

transfer function with a constant gain and a time delay, YpfjoO & Ke ^ 



0 .2 .4 .6 .8 1.0 

Time shift, x, sec 


As noted previously, the value 
of X that will minimize the identi- 
fication error is dependent on the 
effective time delay of the pilot, 
t . For these data, the procedure 
described previously was used to 
determine Tp and, thus, select X. 
With this procedure, Tp was esti- 
mated by fitting the transfer func- 
tion to the plots for several values 
of X. These results are illustrated 
in figure 14 where the estimated Xp 
values are presented as a function 
of X. It is seen that X is equal 
to the estimated time delay, Xp, 
at X ^ 0.6 sec. Therefore, X = 0.6 
sec was selected for use in this 
identification analysis . 


One promising feature in analyz- 
Figure 14.- Comparison of estimated Tp with ing the flight data is that the esti- 

X; example 3. mated describing functions are 

5 The describing function for the jet control system Y c (jw) was 
estimated from the flight data. 
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relatively insensitive to the exact value of the time shift, A. For the 
estimated pilot describing functions (e.g., fig. 13), the plots are approxi- 
mately the same for values of A from about 0.3 to 0.7 sec (the estimated x 
remained the same at about 0.6 sec). It appears that the exact value used for 
A is not critical in this application of the identification technique. 

The estimated describing function for A = 0.6 sec (fig- 13) represents 
a constant gain system with an effective time delay. This result, although 
not directly comparable to the results from previous studies, appears reason- 
able. For instance, with a rate command system, which approximates the control 
system used in this situation, reference 1 has shown that the pilot describing 
function will be essentially a constant gain system with a time delay. The 
value Tp for the three-axis flight data is higher than the value from the 
single-axis data in reference 1. However, other studies such as reference 13 
have also shown higher values of x p when the pilot is involved in the com- 
plete task of monitoring the instrument panel and controlling about three 
separate axes. 

Further analysis of this flight data is presented in appendix B. This 
appendix illustrates how the describing function of the pilot/control combina- 
tion can be identified using the technique outlined in this report. This 
illustration is interesting because it presents the identification of an 
unknown system (i.e., pilot/control system) using only its own internal noise 
source for excitation. 


CONCLUDING REMARKS 


This report has shown that in measuring pilot describing functions the 
identification error due to the correlation of the input signal with the 
pilot output noise can be reduced by shifting the input data during the com- 
puter processing by an amount equivalent to the pilot time delay. 

Both theory and experimental data have shown that the identification 
error can be made small if the autocorrelation function, RnnCO* of the 
internal noise source (pilot remnant) is negligible for x greater than the 
sum of all effective time delays through the control loop (pilot plus con- 
trolled element). This finding has significance in general systems identifi- 
cation because, when these conditions are met, it is possible to measure the 
describing function of a system with feedback using only its own internal 
noise source for excitation. 

Representative data selected from the retrofire portion of the Gemini X 
flight were analyzed using the technique outlined in this report. These 
results demonstrate the feasibility of identifying the pilot describing 
function from routine flight-test records. 


Ames Research Center 

National Aeronautics and Space Administration 

Moffett Field, Calif., 94035, Jan. 15, 1969 
125-19-01-42-00-21 


21 



APPENDIX A 


TIME DOMAIN ANALYSIS AND COMPUTER PROCESSING 


In this appendix, we will first use time domain analysis to outline the 
standard cross -correlation method (refs. 4 and 8). We shall then introduce 
the time shift X and the computer processing equations used for the results 
in this report. The reduction in identification error due to X will then 
be pointed out using these time domain equations. And, finally, we will 
discuss a modification of these computer processing equations to account for 
any data bias. 


Standard Cross-Correlation Method 

Let us consider analysis in the time domain in which the linear input- 
output relationship can be expressed in terms of a convolution integral 

c(t) = J h p (x)e(t - x)dx + n (t) (Al) 

o 

The hp(x) is the pilot impulse response function that is. assumed to be zero 
for t < 0 (i.e., hp(x) is a real system) and also zero for x > t m (i.e., a 
finite memory time, tj^) . A simple discrete approximation of equation (Al) , 
to allow digital computation, is 

M 

c(k) = At hp(m)e(k - m) + n(k) (A2) 

m=o 

where At is the discrete sampling time. The set of equations (A2) can be 
written in vector-matrix form as 


where 


£ = Ehp + n 


e(k 0 ) e(k 0 - 1) . . . e(k 0 - M) 
e(k 0 + 1) 


E = At 


(A3) 


e(K) 


e(K - M) 
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|h ^ 
n P 


c(k Q ) 


n(k 0 ) 


h ™ 
n P 


c(k Q + 1) 


n(k Q + 1) 

lip ~ 

h (2) 
P 

c = 

c(k 0 + 2) 

n = 

n(k 0 + 2) 


h p (M) 


c (K) 


n(K) 


An estimate of hp can then be made, using standard least squares, by the 
formula 

h p = CE t E) _1 E T c (A4) 

T 

We should point out that the matrix to be inverted, E E, contains terms 
that represent discrete measurements of the autocorrelation function Ree( T )> 
and that the vector E^£ contains terms that represent discrete measurements 
of the cross-correlation function Rec( T )- For instance, the vector E^c can 
be written in terms of the cross- correlation function as 



K 




Y eCk)c(k) 

k=k 0 


R ec C°) 


K 




) eCk - l)c(k) 


RecCl) 

_ X 

k=k G 



E c = 

• 


• 


K 

r — i 


• 


\ e(k - M)c(k) 


R ec (M) 


k=k 0 




(AS) 
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Use of a Time Shift in Computer Processing 

The time shift A is introduced into the computer processing by shifting 
the input data e(k) a discrete number of time shifts L, where A = L At. 

The linear pilot model is then expressed as 

M 

c(k) = At ^ hp(L + m)e(k - m - L) + n(k) (A6) 

m=o 

This form assumes that the impulse response hp is zero for time less than 
A. This form also assumes a memory time of t m = A + M At. For the results 
in this report, M = 9 and At = 0.05 sec. (Larger values of M were also 
tried with no significant changes in the results.) 

Using the least-squares formulation, the impulse response function of the 
pilot was determined by the following matrix inversion on a digital computer: 


r 


K 

r — i 

_K 

K 

r — i 

P 1 

_ K 

h p (L) 


) e 2 (k- L) 

) e (k-L)e (k-l-L) 

. . ) e (k-L)e (k-M-L) 


£ e(k-L)c(k) 



k=k 0 

k=k 0 

k=k D 


k=k Q 



K 

K 

K 


K 

hp(L-l) 

= catr 1 

Y e (k- 1-L)e (k-L) 
k=k 0 

V e 2 (k- 1-L) 

U 

k=k 0 

. y e (k-l-L)e (k-M-L) 

k^k 0 


) e (k-l-L)c (k) 
k=k 0 

h p (UM) 


K 

) e (k-M-L)e (k-L) 

K 

) e (k-M-L)e (k-l-L) , , 

/ s 

J 

i 

CM 

<L> 


JC 

) e (k-M-L)c (k) 



_ k =k 0 

k=k Q 

k=k o J 


1 

o 

1 


This time domain solution was further transformed into the frequency domain 
using the following approximation for the Fourier transform: 


M 


YpCjw) = e" L AtJ ’“ At ^ hp(L + m)e 


-m Atjw 


(A8) 


m=o 


-Ajw 


Y m O) 


Reduction in Identification Error With Time Shift 

In order to show, using time domain analysis, that the time shift X 
reduces the identification error, we can write equation (A4) as 


h p = h p + [ETE] _1 E T n 


(A9) 


error 
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The identification error, shown above in vector form, is due to the correla- 
tion of e(k) with n(k). The terms in the vector can be regarded as 

discrete values of the cross- correlation function R en (T). 


K 

r— i i 



) e(k)nCk) 


R en d) 

k=k D 

K 

y e (k - l)n(k) 


R enCD 

kXk D 



K 

\ e(k - M)n(k) 


R en (M) 

_k=k 0 


L _ 


(A10) 


If R en (x) is nonzero for t >_ 0, then the terms R en (m), m >_ 0, will be non- 
zero and there will be an identification error. 


Now, 

equation 


introducing a discrete number of time shifts L, such as used in 
(A7), the vector E T n becomes 


K 

r — i 



/ e (k - L)n(k) 


R en(« 

k=k 0 



K 

r 1 



\ e(k - 1 - L)n(k) 


R e„( L * 1) 

k=k Q 


* 

K 

r — i 



y e(k - M - L)n(k) 


R en C L + M ) 

Jc=k 0 




(All) 


We can see that this use of the time shift removes the terms R en (m), 
0 m < L, and adds the terms R en (m), M < m <_ M + L. The terms that are 
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added generally are smaller than the terms removed; thus, the use of a time 
shift L, or the equivalent X , reduces the identification error. (Time shift- 
ing will not significantly alter the matrix E T E.) Further note that the 
identification error will be zero if R en (m) is zero for values of m >_ L. 


Computer Processing With Data Bias 

The data c(k) and e(k) obtained from routine flight tests will usually 
contain some type of long-term variation about which the short-period dynamics 
are to be estimated. (See, for instance, the flight-test data in fig. 12.) 
There may be a bias or drift in e(k) because of instrumentation errors or 
because the pilot may be controlling about some nonzero value. There may be 
bias in c(k) because of instrumentation errors or because the pilot must hold 
some nonzero offset with the controller (e.g., if the vehicle is out of trim). 
In the analysis of the flight-test records for this report, a bias term and a 
drift term were added to the foregoing formulation so that 

M 

c(k) = b 0 + bjk At + At ^ hp(L + m)e(k - L - m) + n(k) 

m=o 

A As 

The estimated constant bias term b 0 and estimated drift term b x , along 
with the estimated impulse response function hp(m), were determined by the 
least-squares solution. This new formulation required the inversion of an 
M + 3 square matrix instead of the M + 1 square matrix shown previously in 
equation (A7) . 
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APPENDIX B 


IDENTIFICATION OF THE PILOT/CONTROL DESCRIBING FUNCTION 

The describing function of the pilot/control combination has been used 
extensively (e.g., ref. 1) to analyze the closed-loop dynamics of piloted 
systems. This appendix illustrates the use of the technique outlined in this 
report to identify the pilot/control describing function. This illustration 
will use the Gemini data presented previously in figure 12. 

A 

The describing function Y x (ju)) to be identified is shown in relationship 
to the pilot and the control system in figure 15. This describing function 

Pilot 



Figure 15.- The technique for identifying the pilot/control describing function. 

was measured 1 between the attitude error signal e(t] and the attitude rate 
signal iRt) . The Bode plots obtained from these data are presented in 
figure 16. Curves of the estimated describing function Y x (ja)) are shown for 
X = 0 and 0.7 sec. Also shown for comparison is the Bode plot for the nega- 
tive inverse of the feedback path, -jco. The theory predicts that, for X = 0, 
the estimated describing function will identify -jo) because, in this case, 
there are no disturbances external to the measured describing function. As 

l It is interesting to observe that the identification of the 
pilot/control describing function required only a single channel of recorded 
data. The output signal used in the computer processing was the recorded yaw 
rate, ijj(t). The input signal was also determined from this same recorded 
data as 

e(t) = -J KOdT + bias + drift 
o 

where the unknown bias and drift are accounted for by the method shown in 
appendix A. 
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Figure 16.- Identification of pilot/control 
describing function. 


shown in figure 16, the curve for 
X = 0 essentially coincides with the 
curve for -jw. 

For increasing values of X, the 
estimated describing function tends 
away from the curve of - jw. Any 
value of X between 0.4 and 1.0 sec 
resulted in approximately the same 
Bode plot as shown for X = 0.7 sec. 
The estimated describing function can 
be approximated by a describing func- 
tion with a constant gain K x and a 

time delay t x ; Y x (jw) K x £ Tx ^ aj . 

The estimated values for t x are 
presented in figure 17 for several 
values of X. We can see that X is 
equal to the estimated t x at 
X ^0.7 sec, so X = 0.7 sec was 
selected for this analysis. 


From the estimate Y x (ja)), we 
can determine the describing function 
for pilot/controlled element combina- 
tion, YpY c (jcj). The describing func- 



Figure 17.- Comparison of estimated x x with X. 


tion Y x can be combined with the 
integration, 1/jo), as shown in 
figure 15, to obtain 


YpYcOu) 


Y x (ja)) 

jw 


For the results from figure 16, the 
estimated describing function is 

1.3e s0 that 


Y x (jcu) 


YpY c (ju0 


1.3 e ~°‘ 7JaJ 

jio 


This result appears reasonable. From 
previous studies such as reference 1, 
it has been shown that for a variety 


of controlled elements the pilot will control so that 


YpY c 


(jw) 


K e 

X 




jto 


This form is the same as found in the actual flight results. The actual value 
for the gain (a crossover frequency, K x , of 1.3 rad/sec) is lower than pre- 
dicted in reference 1 and the value for the effective time delay (t x & 0.7 sec) 
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is higher than predicted in reference 1. Again, it is reasonable to expect 
(see refs. 10 and 13) that these differences can be attributed to the fact 
that in reference 1, the pilot was controlling only a simple single-axis task, 
whereas, for the actual flight data, the pilot was controlling about three 
axes and monitoring the complete instrument panel. 
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